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Abstract. In this article, we present a new approach to averaging in non-Hamiltonian systems 
with periodic forcing. The results here do not depend on the existence of a small parameter. In 
fact, we show that our averaging method fits into an appropriate nonlinear equivalence problem, 
and that this problem can be solved formally by using the Lie transform framework to linearize 
it. According to this approach, we derive formal coordinate transformations associated with 
both first-order and higher-order averaging, which result in more manageable formulae than the 
classical ones. 

Using these transformations, it is possible to correct the solution of an averaged system by 
recovering the oscillatory components of the original non-averaged system. In this framework, 
the inverse transformations are also defined explicitly by formal series; they allow the estimation 
of appropriate initial data for each higher-order averaged system, respecting the equivalence 
relation. 

Finally, we show how these methods can be used for identifying and computing periodic 
solutions for a very large class of nonlinear systems with time-periodic forcing. We test the 
validity of our approach by analyzing both the first-order and the second-order averaged system 
for a problem in atmospheric chemistry. 
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1. Introduction 



The theory of averaging falls within the more general theory of normal forms; it arises from 
perturbation analysis and is thus formulated usually for equations which contain a small param- 
eter e. A vast literature treats this case (e.g. [21 [TUl CGS [TH [23 E01 1321 US], and references 
therein). In this paper we develop a time-averaging theory for the case of time-periodic nonau- 
tonomous systems without an explicit dependence on a small parameter. The theory is presented 
here formally, subject to a conditional statement to be explained later. It is shown to apply 
numerically to our test problem, which does contain a "large-amplitude" periodic perturbation. 

Among the methods of derivation of averaged systems, those based on Lie transforms offer 
an efficient computational framework (cf. 0131], or more recently [36]). The rigorous use of 
such techniques, however, is confined so far to the construction of near-identity transformations 
expanded in a small parameter. 

In order to overcome the main issues raised by the lack of e-parametrization, we consider the 
problem in terms of "differentiable equivalence" between the time-dependent vector field and a 
corresponding averaged form, allowing us to put the problem in an appropriate Lie transform 
framework. 

In Hamiltonian dynamics, G. Hori |21| was the first to use Lie series for coordinate transfor- 
mations, while A. Deprit [7] introduced the slightly different formalism of Lie transforms; the 
two were shown to be equivalent by J. Henrard and J. Roels |19j . They developed a new, efficient 
way to construct canonical transformations for perturbation theory in celestial mechanics. The 
Poisson bracket in their formulae can be replaced by the Lie bracket to obtain a representation 
of general, non-canonical transformations [5| 116} [23], Nowadays, Lie transforms are applied in 
many research fields, such as artificial satellites [34], particle accelerator design [9] and optical 
devices [12] . to mention just a few. 



The object of the present article is to perform averaging analysis for a T-periodic A^-dimensional 
system (N > 2, T > 0) of ODEs 



where f2 is an open subset of M . The precise framework is given in £J2J 

In the present paper, we consider the averaging process via an appropriate notion of differ- 
entiable equivalence (Definition 12. ip . This equivalence concept allows us to determine the un- 
derlying time-dependent family of diffeomorphisms, as the solution of the nonlinear functional 
equation (|2.6|) . This equation plays a central role in our approach, as it allows one to generalize 
the Lie transform formalism. We put it into a pullback form that falls into this formalism via 
suspension (Proposition 12. ll and Appendix 1). This process allows us to obtain a computable 
formal solution of equation (|2.6|) (Theorem I2.ip . based on the solvability of a recursive family 
of linear PDEs, typically known as Lie's equations [5] 118} 135]. It is shown in §2.21 and §3.21 that 
Lie's equations can be solved easily in our framework. 

This approach is consistent with other classical equivalence or conjugacy problems, in which 
the so-called homological equations linearizing such problems are a cornerstone in proving 
Anosov's theorem or Hartman-Grobman theorem (cf. [2] for details). A recent proof of Kol- 
mogorov's theorem on the conservation of invariant tori also used this classical framework [22] . 
The essential difference between our approach and the classical one is that (|2.6h involves a "con- 
jugacy" (see Definition 12.11 for the right notion) between a nonautonomous field in the original 
problem and an autonomous one in the transformed problem. The formal developments in 
$2] to $4] are contingent upon Lemma 6.1. This Lemma will be proven in a subsequent paper 
(in preparation), and permits a rigorous proof of the existence of a solution to equation (|2.6p . 



(1.1) 
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The numerical results in [j5] and £[6] demonstrate the usefulness of the present approach and the 
plausibility of the forthcoming rigorous results. 

Our approach provides, furthermore, simpler formulae for higher-order averaged systems, 
even in the classical context of perturbations scaled by a small parameter e. This is shown in 
Proposition 13,11 and is due to the suspension associated with solving (|2.6|) : see ^3.2^ §3.31 and 

PZ3I31. 

We show that the classical change of coordinates on the initial data respects the differentiable 
equivalence relation (|2,6p at time t = 0. It arises here in a more natural way from the generator 
of the family of inverse diffeomorphims at time t = 0; see Proposition 14,21 and |17j . 

We demonstrate that the Lie transform framework permits gaining CPU time with respect 
to the standard numerical methods by, loosely speaking, integrating the averaged system with 
a larger step size; this is followed by performing efficient corrections with a step size adapted 
to the oscillations of the forcing terms. We test the validity of this numerical approach by 
computing the second-order terms in our atmospheric-chemistry problem. More precisely, we 
build a second-order averaged system in this test problem and the relevant correction in the 
same framework, after applying the averaging-correction method at the first order. 

In £}2]we describe how Lie transforms allow us to resolve formally a differentiable equivalence 
problem between a time-periodic system and an autonomous one. In we give a formal 
algorithm of computing higher-order averaged systems based on Lie transforms and related 
corrections; furthermore we compare our results to the classical ones. In §H we clarify how to 
compute initial data which satisfy the differentiable equivalence relation. Numerical results and 
advantages of this technique are discussed in ^5] for a simple model of atmospheric chemistry. In 
£j6j we comment on the use of these methods for constructing periodic orbits, in nonautonomous 
dissipative systems with periodic forcing. In the Appendix 1, the Lie transform framework is 
presented as a tool for solving pullback problems in general. 

2. The Lie transform setting and the differentiable equivalence problem 

2.1. Fitting the differentiable equivalence problem into the Lie transform formalism. 

For all the rest of the paper fi will be an open subset of Mr, N will be a positive integer, and 
T will be a positive real. 

For subsequent computations we shall consider Q fc (f2) as the class of continuous functions 
/ : M + x — > W N which are T-periodic in t for each x € f2 and which have /c-continuous 
partial derivatives with respect to x for x £ and t in M + . More precisely we will work with 
Q°°(f2) = Hfc>o Q k (^) wnen we deal with Lie transforms as often it is requires in this framework 
for dependence on the spatial variables (cf. Appendix 1). 

There are known results on the boundedness and the global existence in time of solutions of 
ODEs. For instance, by an application of Gronwall's lemma it is easy to show that if there exists 
a,P G C (M + ,R + )nL 1 (IR+,lR) such that \\Y(t,x)\\ < a(t)\\x\\+ f3(t), for all (t,x) £l+x(l, then 
every solution of x = Y(t, x) is bounded. Thus, using the continuation theorem (cf. theorem 2.1 
of p3] for instance), we can prove that the solutions of the above equations are global in time. 

Taking into consideration this fact, we make the following assumption. 

(A) All solutions remain in Q, for all considered non-autonomous and autonomous vector 
fields defined on CI. 

Henceforth, for all 1 < k < oo, we introduce V k {Cl) := Q k (Cl) n {/ : R + x CI Cl} and 
C k (Cl) := C k (Cl,Cl). 

We describe some concepts related to the subject of this paper, dealing with differentiable 
equivalence between a nonautonomous vector field and an autonomous one. 
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Definition 2.1. Let r be a positive integer, let Y £ V r {Vt) and Z £ C r (£l), then Y and Z are 
said to be 'P^ft- equivalent (k < r), if there exists a map $ E V k {£i) such that for all t € R + , 

<3?j := <3?(i, •) : SI — > $7 is a C k -diffeomorphism, which carries the solutions x(t,xo) of x = Y(t,x), 
for x(0,xq) = xq varying in 0, into the solutions z(t,zo) of z = Z(z), in the sense that: 

(2.1) x(t,x ) = Mz(t,% 1 (xo))), for allteR + ,(z = % 1 (x )). 



Using the notations in this definition, we shall say that if (|2.ip is satisfied for a pair (x, z) 
of solutions, then x and z are in V\^t -correspondence by <]>, or, sometimes that they are in 
V^if f -correspondence by (<E> t ) teR +. 

Furthermore, we will often use for the underlying transformations the functional space 

V k {Q) := {$ € V k (n) such that, for all felR + ,$ ( :fl^!] 

(2.2) is a C fc -diffeomorphism, and t -)• ^j -1 is C 1 }, 

where the smoothness assumption on the map t — > Q^ 1 , will be apparent from the proof of 
Lemma 12.11 

Remark 2.1. It is important to note here that the V^t j-- equivalence of Definition \2.1\ does 
not conserve the period of the trajectories although it preserves parametrization by time and 
sense, and, in this meaning, realizes a compromise between the classical notions of conjugacy 
and equivalence (e.g. [13 \). This will be essential in £J21 

Let us recall the classical definition of pullback of autonomous vector fields. 

Definition 2.2. Let U and O be open subsets ofM. N . Let X be a vector field of class at least C° 
on O, and let $ 6 C l {U, O) be a C 1 -diffeomorphism. The pullback of X on O by <3? is the vector 
field defined on U as the map 

(2.3) y ($- x * X)(y) := (D$ _1 )($(y)) • X(Q(y)), for all y € U = ^{O). 



Remark 2.2. Sometimes we will employ the notation (<3?*X)(y) := (D<J>(y)) 1 • X(Q(y)) for 
^-UX^y). 

The problem of "P^j-equivalence takes the form of a family of "pulled back" problems via 
the 

Lemma 2.1. Let Y and Z as in Definition \2.1\ Let x : M + — >• Q. a solution of x = Y(t,x) := 
Yt(x) ,x(0) = xo 6 SI, and z an integral curve of Z. Then, x and z are in V^f f- correspondence 
by $ € V^{VL) if and only if z is solution of the IVP: y = * Y t )(y) + dt&^i^t (y)) , y(0) = 

Proof: The proof of this lemma is an obvious application of the chain rule formula and Definition 
EU □ 

Let us now introduce the projection 

(2.4) 7T : 



(t, x) — > x 

We define also a map, defined for each t G R + , as 
(2 ' 5) lt -\x^(t,x) 
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which will be T-periodic in time t. 

According to Lemma 12.11 if there exists a map <3? € which satisfies the following 

nonlinear functional equation, (called in this work equivalence problem) 

(2.6) * Y t + d^ 1 o $ t = Z, for all i G R+, 

then every integral curve of Z is in "P^j-correspondence by <3? to a unique solution of: x = 
Y(t,x), that is the vector fields Y and Z are V\* ^-equivalent. 

The problem is then to solve in $ € V%(Cl), the equation (|2.6p for a given pair of fixed vector 
fields Y € V r (Q) (non-autonomous) and Z 6 C r (S7) (autonomous). 

In order to show how Lie transforms formalism (e.g., Qj|] and Appendix 1) can be used 
for solving formally equation ()2.6p with r = k = oo, we put this equation in an appropriate 
pullback form by the Proposition 12.11 below. 

For that we introduce some notations. We denote by Y the so-called t-suspended vector 
field associated with (jl.ip written in terms of the enlarged vector of dynamical variables, x = 
(t,x\, ■ ■ ■ , xn) T , thus if x is a solution of x = Y(t, x), then x satisfies 

(2.7) ^ = Y(x) where Y = [1 R , Y] T . 

Note that by assumption (A), every solution x of (|2.7|) is contained in 17 := {(t,x) € M + x $7}, 
for the rest of the paper. 

Remark 2.3. iVote that for the rest of the paper the t-suspended form of a (time- dependent 
or not) vector field X on will be defined as X = [l^,X] T , and the flat form of X will be 
defined as [0^,X] T . The same notation X will be used for the latter, by abusing the notation. 
No confusion will be made with regard to the context. 

Proposition 2.1. Let Z € C°°(f2) be an autonomous vector field and let Y E V°°(£l), with Z 
and Y as their respective t-suspended forms. Let 6 V^{VL). If satisfies Q^ 1 *Y = Z , for 
all t > 0, then the family of C°° -diffeomorphisms (it o Q t ol ( ) (€ |+ ; acting on ft, satisfies the 
equation \2.b}) . The converse is also true. 

Proof. The proof relies on basic calculus and the definitions listed above. □ 

We illustrate now the strategy used to solve equation (|2.6p for Z = Y := [l,^] with Y 

being the classical first averaged system, namely ^ Jq Yds. The same reasoning will be used 
for higher-order averaged systems in §3j By Proposition I2.lt the- problem of finding a solution 

of (12. 6p falls into the Lie transform formalism if Y and Y are given by r-series (cf. Appendix 1 
and Theorem 7.1). This point of view, according to the solvability of Lie's equations, will then 
allow us to solve formally ()2.6p in §2.31 

Starting from that point, denote Y by Y ave . In order to write Y ave and Y as r-series we 
introduce a natural auxiliary real parameter r, by the following expressions 



(2.8) Y aveiT := [1 M , R ^] T + r • [0 R , Y] , 

(2.9) % := [1 R , R ^] T + r • [Or, Yf = Y (0) +r-Y 1 

with the obvious definitions for the vector fields Yq ^ and Y^ , and 

(2.10) y fc (0) = R iv+i, for all integers k>2. 



6 



M.D. CHEKROUN, M. GHIL, J.ROUX AND F. VARADI 



Consider (p T i the semiflow at r generated by a r-suspended smooth field Gt for each non- 
negative t, i.e. tp T) t is the solution at r of 

(2-H) S G t {i = GS = G(U) ,£ = (r,0 € R^ 1 , 

with a formal expansion of G(i, £) in powers of r given by 

(2.12) = l,G(t,0 =J3- r G n (t,0 J 



n>0 

where 

(2.13) G (t, = [1, G (i ; 0] T and for n € 1? + , G n (t, £) = [0, G n (t, £)] T , 

with G and G n smooth vector fields on R . 

Remark 2.4. T/ie generator system would be denoted by another notation, the r-suspension 
and the t-suspension being not equivalent generally. We slightly abuse the notation. 

Applying the discussion of the Appendix 1 with p = N + 1, A T = Y T and B T = Y ave ^ T for 
all t, and utilizing the fact that Lie's equations are independent of r, we obtain that if there 
exists (p T) t such that (ip T j)*Y T = Y avetT for every non-negative real t, then by Corollary 7.2 with 
H Tjt = G Tt t, Lie's equations 

f y (0) = [1m, (M T , 
(2-14) y (1) = [o R ,F] T , 

[ Y^ m) = [OV+i] for all m G Z + \{0, 1}, 

are solvable for the G n j- 

Note that Y (resp. Y ave ) corresponds to r = 1 in (j2.9|) (resp. (|2.8p ). so if (ipij)*Y = Y ave for 
every non-negative real t, then by Proposition 12.11 the one-parameter family (ir o ipi >t o I t ) t€ ^+ 
solves the equation ()2.6p . We solve Lie's equations (|2.14|) in the next section and we will show 
how this solvability will allow us to compute formally a solution of (|2.6j) in §2.3. 

2.2. Solving Lie's equations associated to the equivalence problem (12. 6|) - first-order 
averaged system. We treat in this section the problem of solvability of Lie's equations asso- 
ciated with (|2.6p for Z = Y ave . 

In that respect, let M be the map from x to the set of the vector fields of M Ar+1 defined 

by 



(2.15) 



M(j + k,j) =Y { k j \ \/(k,j) G 
M(i,j) = r n+i i<j, 



where the terms Y{, are calculated from those of (|2.9p and (|2.10p by the recursive formula (|7.5p . 
We have the following result: 

Proposition 2.2. For every integer j greater than or equal to two and for every integer I 
belonging to {0, .., j — 2}, we have M(j,j — I) = 0. 

Proof. We prove this fact by recurrence. First we note that M(2, 2) = Yq = by ()2.14p . If 
j = 3, then the recurrence formula (|7.5|) applied for n = and i = 2 for the vector field Y and 
re-written for M, gives M(3,2) = M(3,3) - L go M(2,2) but M(3,3) = M(2,2) = because of 
(HHD, so M(3,2) = 0. 
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We make the following hypothesis, V(j), until a j > 3, 

V(j) : M(j - 1, (j - 1) - I) = 0, for all / G {0, .., j - 3}. 

Then, as a consequence of this predicate, we get 

Vp E {0, .., /}; M((j -l)-Z + p,0-l)-0 = 0. 

We show the heredity of the predicate. The formula (|7.5p for n = I and i = j — (Z + 1), takes 
the form 

p=0 

so is re-written, for M, 

M(j, j - = M(j,j - (I + 1)) + J2 C?L 5[ _ p M(j -{l + i) +P) j-(l + 1)), 

p=0 

and, because of the predicate, 

(2.16) M(j,j -l)= M(jJ -I- I) for all I G {0, .., j - 3}. 

But for I = 0, M(j,j) = Yq^ = according to (|2.14[) . so M(j,j — I) = for every integer I 
between and j — 3. Making I = j — 3 in ()2. 16j) we get M(j, 2) = and therefore V(j + 1). □ 

We make a useful remark for the clarity of some coming computations. 

Remark 2.5. Let Y be the vector field associated with the system (1.1), with an expanded t- 
suspended form given by (2.9) and (2.10). Let G be a vector field on M. N with G given on the 
model of (2.12) and (2.13). Let Y^ 1 be the terms calculated by the recursive formula (7.5) with 
A = Y and H = G. Then ir o yj™ 1 ^ = Yj m ^ for every (m,j) G r L + x Z+. 

This last proposition allows us to solve Lie's equations easily in our case, which we can 
summarize in the following 

Proposition 2.3. Suppose that (pn generated by Gt given by \2.12}) . satisfies {ip\t)*Y = Y ave 



for all t > 0. Then the sequence (Gj,Yj = a^ 1 ' ,Yj L> )jezi (with a^ 1 G M.) is entirely 

determined by the sequence (Yj^)j&* + given by h2.10\) and h2.9\) . 
More exactly we have, for every positive integer j, 

(2.17) af = 0, Gj(t) = J (JLa^Y - Y^)ds , Vt G R + 

modulo a constant vector, and, 

k=j-l 

( 2 - 18 ) ^ = - E c U l gy£U- 

k=0 

The initial term (GqjYq 1 ^) is defined by the first Lie equation and \2.1J$ . 

Proof. First we prove that Go (and Go) is determined as follows. 

By definition of the Lie derivative in R , the first Lie equation (cf. Appendix 1) takes the 
form 

(2.19) y (1) = L 5 Y^ + if> = DY^.G - DG .Y® + 



T 



.,(!) 
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where Go is the unknown. Since DY (0) = [0 jj(jv+i)x(jv+i)], (I2.19D becomes 

(2.20) y (1) - Y} 0) = -DG .Y^ 0) . 
where 

(2.21) DG .Y^ 0) = [0 R ,d t G ] T . 

Hence, the equation (|2.19p for Go becomes, taking into account the expressions of y} 0) and 
y (1) given by ([2JjD and (gZHD , 

(2.22) d t G (t,0 = Y(t,0-Y{0- 
which gives Go by integration. 

We obtain as follows. By (|7.5p we get 

9(2) _ ^(l) , _ ^(1) 

but y (1) is given m (I2TT41) and y o (2) is equal to zero. Thus we get y} 1] = -Lq F (1) , which yields 
(gJSD for j = 1. 

Computing Gi by (|7.5|) gives: 

and, after easy manipulations as Y^ = K iv+i by (I2,10p . 

= 0, 

d t G 1 = L Go {Y + Y), 

yielding (|2.17p for j = 1 after integration. 

We consider the following predicate, H(j), for j > 1: 



af ] = 0, Gj(t) = fiVLe^Y - Y} l >)ds , Vt G M+, 
W(j) : <( and 

y(l) — _ \^k=j-l r k j _ y-(l) 

We prove %(j) by recurrence. The preceding computations show H(l). Let j > 2. Suppose 
H(j - 1). Again (|7.5p gives, 

(2.23) 5f = y w + L 5 /r + £ ^_/ fe (0) 

fc=l 

and, as previously, we have 

(2.24) L d Y^ 0) = [0 R ,-d t G j f . 

The terms Y" fe {0) are equal to zero for all k > 2 (see (gZED ), so, using (2.12), (2.13), (g^gj) , and 
y^ -* = y ; we get by integration the two first assertions of noting that the first component 

of the r.h.s of (|2.23p is equal to zero (that gives aij = 0). 

Finally, we obtain the last assertion of H(j) as follows. From (|7.5p . we have for every positive 
integer j: 

(2-25) Y} 2 \ = Jf > + £ C^L d Y^_ k , 

k=0 
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but by Proposition 12.21 we get in particular that M(2 + j — 1, 2) = = for every positive 
integer j, so we deduce the last assertion of ~H(j) from (|2.25p . Thus rl(j) is a necessary condition 
of T~L(j — 1), that completes the demonstration. □ 

2.3. Formal solution of the equivalence problem (12.61) . We describe in this section a 
procedure for obtaining formal series representation of a solution of (|2.6p . assuming that Lie's 
equations are solved. In other words, we assume, for this section, that Lie's equations associated 
with the problem of equivalence under consideration are solved. 

Introduce for any sufficiently smooth vector fields F, H : (r, £) € R JV+1 — > W N , 



(2.26) 



dF 

A H F = — + D^F-H, 



where D^F stands for the usual Jacobian of F in the local coordinates £ € 
Naturally, represents the iterated operation 



A™ 



o A 



H 



The following lemma is the first step toward efficiently describing "PJ^-j-correspondence be- 
tween a solution x of (11. ip . and one z of y = Z(y). 



Lemma 2.2. Let 4>ij generated by the r-suspended vector field Ht 
Z for all t > 0, then tt o <fii t o Tt is given as follows 



{l,H t ) T such that (<fii,t)*Y 



(2.27) 



7T O . 



M ol t = Id RN +H , t + A Ht H t\ T=0 ' f° r al1 1 € 

i>l ^ >' 



Proof. This lemma is a consequence of the Taylor formula. Let z be an integral curve of y = Z(y), 
and x the solution of (jl.ip being in P^^-correspondence with z by (tt o cfi 1 j oX t ) ieK +. We have 
formally for r = 1: 

(2.28) (</> M o Xt)(z(t)) = 0(1, t, (t, z(t))) = (t, z(t)f + d 4- ^ ' r/ "~ J ' 



dr 



C n> 



' (n + 1)! dr n+1 



C 



where C means the triplet (r = 0, t, z(t)). 

Consider the projection tt on the last components (projection which commutes with the 
operator 4-). If we compose (I2.28P to the left by tt we obtain: 



(2.29) 



(vr o 1>t oZ t )(z(t)) = (tto 4> ltt )(t, z(t)) 

d(rr o 4>) 



z(t) + 



dr 



d n+l (TT 



^ (n+1)! dr n+1 
n>l 



taking into account that -£p (tt o 0) | = 
integer n greater or equal than 1 that: 



Ho(t, z(t)), it is therefore sufficient to prove for every 



(2.30) 



d n+1 (lT- 



dr n+1 



C 



A n Ht H t (r,0\ c . 



We prove (|2.3U|) . Using the chain rule it is merely an exercise in basic calculus. 
By assumptions, 



(2.31) 



H, 
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Let us introduce the obvious notation <fi T j = (Idjt; (p^ t ; ...; 4>^ t ) T and, denoting the i th com- 
ponent of Ht, T by Hf T , the notation Ht :T = (l,Hj T , ...,H[ f T ) T . 

Moreover, let us recall that the classical pullback by a diffeomorphism $ of a map / : M. N — > R 
is defined by: 

(2.32) **/ = /o*, 
thus (|2.3ip can be viewed component-wise as: 

(2.33) = (<M*#lr, for all i G {l,...,iV}. 

CLT 

Noting that (<fi T j)* H% T (-) = Hf T o <ft Ti t(-) = H\{t, <Ar,t( - ))> by application of the chain rule, 
taking into account (|2.32p . we obtain finally the key formula: 

(2-34) jj.((^ >t )*flJ T ) = (cp T , t y(d r Hl T + (VHl T ,H t>T )), 

where V stands for the gradient and (•, •) for the usual Euclidean scalar product on R N . 
By introducing the following operator, which acts on families of maps from R^ to R: 

(2-35) T Ht (-) = d T (-) + {V(-),H t ), 

we can write ()2.34p in the form: 

(2-36) A ((<M *^ t ) = (^, t )*T^ )T (flJ T ), 

and thus, by induction it is clear that, for every integer n, 
(2-37) _( ((M *tf* T ) = (0 T)t )*T^ r (H t V), 

which, according to (12. 33p . leads to 

(2-38) ^r«t) = (W*T&«, T (^r), 

and taking this expression at r = we get, because 0o,t = Id^N+i: 

(2-39) ^+i(€,t) \r=o= T% T (Hlr) \r=o, for alH G {1, • • • , N}. 

Let F be a map from R N+1 to R N . By definitions (f!Q2"j) and (i236|) we see that 
(2.40) A Ht F=(T Ht F\... ,T Ht F N f, 

so if we apply (12.40P to F = H t , taking into account (|2.39p . we obtain: 

(2-41) d^ {7T ° ^ |t=0= A ^ Ht |r=0 ' 

which is (|230l) . 

□ 

Consider now and time-dependent vector fields acting on R^, both admitting r-series 
representations. Omitting the dependence on t for the vector fields in the expansion of Ft and 
Ht, we consider the following recurrence formula: 

n 

(2-42) i^+H = + £ CnDiFf ■ H n ^ k ; for all i € Z + , 

fc=0 
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initialized by Fn = F U)t , F n j being the n th term of the r-series of F t ; and where Fl n -k is the 
(n — k) th term in the T-series of Ht- 
With these notations, we formulate, 

Definition 2.8. Let t € M + fixed. The Hftransformation evaluated at t, denoted by 7~H t ( T )> 
of a time- dependent vector field Ft acting on M> N , by a time- dependent vector field Ht acting on 
the same space, both admitting formal r-series, is given by: 

^ T n+1 . 

(2-43) T Ht (r) -F t = J2 T^Tq! F oS > 



n>0 



where each -Fq™ is calculated using the recurrence formula \2.^2^ , initialized by F\ l t = F n j; F U) t 
being the n th term of the formal r-series of F t . When F = H, the term F^ will be called the 
(n + l) th - corrector. 

Remark 2.6. It is important to note here that 7~H t (0) ■ Ft = 0, allowing near-identity transfor- 
mations for small r (cf. Theorem \2. 



We then have the following theorem: 

Theorem 2.1. If <j) T t is generated by a time- dependent vector field of the form (l,Ht) T on a 
region M. + x 0, for which Ht has a formal r-series, Yl n >o nT-^M' on ^ s re 9i° n f or a ^ T £ I> 
I being an interval of reals such that £ /, then: 

(2.44) vr o (f) T;t o%t = Id RN + T Ht (r) • H u on I x R+ x Q. 

Proof. First of all, note that the derivations in the proof of Lemma [2. 21 show that, for all t € M + 
and t £ I, 

T i+l 

(2-45) vr o <rV,t o l t = Id m + r • H Qjt + — — A Ht H t \ T=Q . 

i>l 

Omitting the time-dependence, simple calculus shows 

n n 

(2.46) A H H = ^(fln+l + °n D t H k ■ H n -k), 

n>0 ' k=0 

and defining the vector fields H^ as 

n 

(2.47) H W = H n+1 + ^ C^D ( H k ■ H n _ k , 

k=0 

following a procedure similar to the one of |19j . but for the operator Ah, we get the recurrence 
formula: 

n 

(2.48) H^ =H® +1 + Y CnD^Hf ■ H n _ k ; for all i e Z+, 

k=0 

where H$ is the n th term of the series 



T 

(2.49) A i HH = J2^H®, 

n>0 

and the initial terms are given by Hn = H n for every non-negative integer n. 
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We conclude therefore that: 




□ 



In practice Theorem 12,11 and Corollary 7.2, permit to determine formally a solution of (j2.6|) . 
For example, when the "target" autonomous field is Y, as shows Proposition 12.31 Lie's equations 
(|2,14j) are solvable, thus (Id^N + 7c? t (l) • Gt)tes+ can De computed using computer algebra 
implementation O [25l [6] and thereby gives a formal solution of (|2.6p (with Z = Y ) , where the 
generator G is determined by (12. 17ft in Proposition 12.31 



Higher-order averaged systems are linked to the way change of variables are performed. In 
the e-dependent case, the usual formulae are obtained by using formulae for arbitrary-order 
derivatives of compositions of differentiable vector functions [11], on the one hand, and by using 
the implicit function theorem, on the other |1Q |, I3Q| . [32] . This form for higher-order averaging is 
not the only one, as it is pointed out in |29[ pp. 36-37]. We describe here other higher-order 
averaged systems both for e-dependent and independent cases, via Lie transforms formalism, 
but different from those of [311 [36] . and given by simple explicit formulae. 

3.1. Formal algorithm of higher-order averaging based on Lie transforms. We de- 
scribe here a way of computing formally the higher-order averaged systems according to the Lie 
transform of the vector field Y S V 00 ^) in consideration. 

In that respect we consider first the parametric standard form x = rY(t, x), where r is a real 
parameter. Then we have the following proposition. 

i n \ 

Proposition 3.1. Let n be an integer greater than or equal to two. Let Y £ C°°(Q) be the 
n th term of the finite sequence satisfying 

(a T ): = Y ^ + T* +1 [/j + i, where t/j+i is an autonomous field in C°°(f2) for all i € 

{0, • • • ,n-l}, with Y {0) =0,U 1 = Y and U = 0. 

Assume that there exists <p Tt t generated by W% = [l,Wt) T having a r-series represen- 
tation such that 

(b T ): For all i E {1, • • • , n — 1}, the i th term Wi of the r-series representation of Wt is 
T-periodic. 

Lf 7r o <f> T>t o Xt satisfies for all t € M + 



3. New higher-order averaged systems and related corrections 



(3.1) 



(vr o T)t oXt) * (rY t ) +dt(iro 4> Ttt o X t ) 



(-7T O (j) Tjt oX t )=Y 



then 



(3.2) 




m=2 



where for m 6 {2, • • • 



n} 




.T m—2 k=l 



(m-l-l) 



(3.3) 




) + C m _ 1 L H / m _ 2 y)ds; 



1=0 k=0 



and Wo = Go, with Go given from h2.22\) . 
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Proof. Let an integer n > 2. If there exists (j) T> t generated by Wt = [1r, Wt] T having a r-series 
representation on the model of (|2.12p and (|2,13p . such that for every positive real t, the map 
7ro0 r t oI ( satisfies (|3.ip . then according to Proposition 12.11 the equivalence problem f|3. If) takes 

the form of the pullback one * ( T Y) = Y [n \ that gives by Theorem 7.1 

L(W t )(r)-(rY)=Y in \ 

i n \ 

where Y is given by assumption (a T ). 



If we applied (7.4) in this case, we get: 

m=oo m 

V —Y (m) - V r m U 



m=oo „, m=n 

ml 

m=0 m=U 



according to the assumption (a T ), where for rn > 1, U m (resp. Yq™ ) is the flat form of C/ m 

(resp. ^ot )' an< ^ * ne ^-suspended form for m = 0. 

By projecting onto the iV-last components of R JV+1 , by taking into account the remark before 
Proposition 2.6 we get for every m € {2, • • • , n}, 

(3.4) ^ m = J_y M 

ml 

giving thereby (3.2). 

(3.5) y (m) =Y} m - 1) + L Wo Y^ m - 1 \ 
expressing Y^ m ^ by (17. 5p we get: 

(3.6) Yf^ = Y^ m - 2) + £ ClL Wl _ k Yt m - 2 \ 

k=0 

and continuing the procedure to an arbitrary rank I G {2, ■ ■ • , m — 1}: 

(3.7) y/^> = y^r 1 " + E c?^^™- 1 - , 

fc=0 

thus collecting, in ()3.5|) . from / = 1 until Z = m — 1 successively the terms Y^ 1 ^ we get: 

m—l k=l 

(3.8) y H = y w + 53 ( E Q^yf- 1 ^). 

Z=0 A:=0 

Noting that y/ ^ = 7r o is equal to K (see (|2.9|) ) and taking into account (|2.10p . we get 

( 3 -9) E C m-l L W m -l-k Y iP = C ln-l L W m -2 Y + L W m -i Y o°\ 

k=0 

and noting moreover that, as in (I2.24p . 

r(0) 



Express the term Yq 71 ' . By (|7.5|) we have for m E {2, • • • , n}: 



m— 1 



(3.10) £w in _ 1 i r o w = ~W 



m—l) 
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the formula (|3.8p takes the form 

m— 2 k=l 

(3.11) y (m) = Y J (Y, c i L ^ Y t^ l) ) + c 1 m - 1 L Wm _ 2 Y-d t w m ^. 

1=0 k=0 

Integrating (13. lip with respect to time, from to the minimal period T, we obtain according 

to the fact that Yq 7 ^ is autonomous by assumption (a T ) and the fact that W m —\ is T-periodic 
by assumption (6 T ), formula (j3.3|) taking into account f)3.4|) . 

Finally, noting that the first Lie equation is the same at each-order (cf. (|3.18p and related 
discussion), we easily conclude that Wq = Gq, where Go is given by (|2.22p . which completes the 
proof of this proposition. □ 

Note that when we talk about the n th averaged system, in any cases, we will consider naturally 
the following system of differential equations: 

(3.12) x = Y {n \x), 
where Y will be given by (I3.2p in Proposition 13.11 

3.2. Correction of higher-order averaged systems and solvability of related Lie's 
equations. We adopt here the procedure described in §2.3l for correcting higher-order averaged 
systems. 

We consider the system of differential equations (jl.ip . and for a given p > 2 a higher-order 
averaged form given by x = Y (x) where Y is given by (|3.2p in Proposition 13.11 for r = 1. 

Consider the diffeomorphism <j>n generated by Wt = [l,Wj] T in Proposition 13.11 Then, 
according to Theorem 12.11 the transformation ir o (f>xj o I t is given formally as the following 
formal series representation 

(3.13) Id MN +T Wt (l)-W t , 

that allows computer algebra implementation [6], and where Wt depends on the averaged system 
f„\ 

Y KF > (cf. again Proposition 13 . X [> . 

It is important to note here that for two given averaged systems, Y and Y , (p ^ q), 
the corresponding generators, say W and V for fixing the ideas, related to each system are 
not a priori equal. Therefore, considering Theorem 2.9 and Definition 2.8, the Wt- and Vt- 
transformation given by analogous formulae to (|3.13p . are also not equal. 

For instance for p = 1 (resp. p = 2), i.e. when W = G (resp. W = T), G (resp. T) being 

defined as the generator of time-dependent diffeomorphisms which transform Y into = Y 

(2) (2) 

(resp. Y = Y + Yq ); the second corrector in 7g*(1) ■ Gt and in 7r t (l) • Tt are respectively 

(3.14) G [ ^ = G ht + DG ,t.G ,t, 
which is given by (|2.42p with F = H = G, i = 0, n = and, 

(3.15) r$ = r 1 , t + Dr 0j t.ro,t, 

which is given by (|2.42p with F = H = T, i = and n = 0. 

For p = 2, we have Tq^ = Gqj (see (|3. 18|) below) and, by applying (|7.5p for i = 1, n = 0, 
A = Y, and H = T = [Ir, T] t , with some easy computations we get 

(3.16) r lit = (D(Y + Y).r -DT .(Y + Y)-Yi 2) )ds, 
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modulo a constant, where Yq is given by (13, 3D with m = 2. 

Same algebraic procedure for obtaining T\ t t and related computations, in the case p = 1, leads 

to 

(3.17) g m = r(z)(y + F).Go-z?G .(y + F))(fc, 

•/ o 

modulo a constant. 

Therefore, we see that the Gt-transformation and the rj-transformation are different. Indeed 
by substituting (3.17) and (3.16) respectively in (3.14) and (3.15), we note that the second 
correctors Gq| and F^~\ are distinct. Same remarks hold for higher-order correctors related to 
two distinct higher-order averaged systems, allowing us to conclude the non-invariance of every 
higher-order corrector. 

We consider now the invariance of the first corrector. Indeed, the relation (|7.5|) for i = 0, 
n = 0, A = Y and Hq = Wq or Hq = Vq yields for the same Lie's equation 

(3.18) d t (7roH 0it ) = Y t -Y 

after projection. This demonstrates that the first corrector term W j is independent from the 
averaged system considered. 

Note that the non-invariance (resp. invariance) of higher-order (resp. first-order) corrector 

Wq U } (n > 1) (resp. Wq ]), was mentioned in [291 P- 36] in the context of general averaging (i.e. 
without the use of Lie transforms). 

We conclude this section concerning the solvability of Lie's equations related to higher-order 
averaged systems. 

Fix p > 2. For computer implementation of (|3.13p . we have to initialize f|2.42|> with H = 
W = F; in other words we have to solve Lie's equations associated with (13. ip . Those equations 
show a structure similar to that studied in §2.21 We outline how to solve them for p > 2, and we 
give the precise results for the case p = 2. Here, the key point of the Lie's equations solvability 
consists in extending Proposition 2.5 associated to the first-order equivalence problem (i.e. for 
n = 1 in (3.1)) to one associated to a higher-order equivalence problem (3.1). Indeed, based on 
(|7.5p . by introducing a map J\f on the model of (12.151) with elements of the diagonal given by 
(3.3), we can show that a proposition analogous to Proposition 12.21 holds, by shifting of (p — 1) 
ranges the first column of vanishing terms in the infinite matricial representation of N . This 
allows us to compute the Wj (i.e. the j -term of the r-series representation of W in Proposition 
3.1), by adapting computations giving the Gj in Proposition 12.31 We illustrate this procedure 
in the case p = 2 for the convenience of the reader. 

In this latter case, we can indeed show the two propositions below. First of all note that, by 
applying (3.3) in Proposition 3.1 with m = 2, we have 

(3.19) y (2) = ~.J (D(Y + F)T - DT .(Y + Y))ds, 

where Tq is given by integration of the r.h.s of (3.18) because of the invariance of the first 
corrector. 

The analogous of Proposition 2.5 associated to the second-order equivalence problem (i.e. 
n = 2 in (3.1)) is 

Proposition 3.2. Let M be a map introduced on the model of 112.15)) , such that M(l, 1) = Yq 1 ^ 
is given by (2.14), Af(2,2) = Y^ 2) = [0,y o (2) ] T where Y Q (2) is given by (8J$) , andN(j,j) = Y (j) 
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are null for every integer j greater than or equal to three. Then for every such integer j, and 
for every integer I belonging to {0, j — 3}, we have Af(j,j — I) = 0. 

Following the model of the proof of Proposition 2.6 that was based on Proposition 2.5, we 
can show, taking into account Proposition 3.2, 

Proposition 3.3. Suppose that ipij generated by Tt = [lR,rj] T satisfies ip^l * Y = Y^ for all 

t>0. Then the sequence (T 'j ,Y^ 1 \y^ 2 ' ') w entirely determined by the sequence (Y^)j^j* 
given by (2.10]) and (2.9). More exactly we have for every positive integer j, 

ft 

(3.20) Tj(t) = / (jLr^Y - Y^)ds , Vt € R+, 

J 

modulo a constant vector, 

(3.21) if = if - £ C^L??} 



(1) _ ^(2) \- n k T „ v(V 

1-fc' 

fc=0 



and 



k=j-l 

(3-22) Y^ C"-iL f Y^_ k , 

k=0 

with every first component ofYj 1 ^ andY^ being equal to zero. The initial term {Tq, Yq , Yq 2 ^) 
is defined by the first Lie equation, \2.1J$ , and (3.19). 

This concludes the solvability of Lie's equations related to the second-order averaged system. 
Such propositions can be extended for higher-order averaged systems, that rules, by induction, 
the solvability of related Lie's equations. 

3.3. A heuristic comparison with the classical approach. We make in this section a few 
comments on previous works and ours as well, and discuss a rigorous setting concerning formulae 
derived in the preceding sections, in the case when r is a small parameter, usually denoted by e. 

As explained in §3.21 Lie's equations are solvable to each order of the averaged systems, 
giving thereby the time- varying change of coordinates by Theorem 12.11 For instance, for n = 2 
in Proposition I3.1| we can solve Lie's equations related to the second-averaged system (cf. 
Proposition 3.3). Consider the system (1.1), if we make the change of variables 

2 

(3-23) x = y + eWot ( y ) + L W lj( y ) j 

(that is a truncation up to order 2 of 7vk ( (1) • Wt), then following the procedure of [HI PP- 
168-169] under sufficient smoothness assumptions, we obtain that y satisfies 

2 

(3.24) y = e Y + € -Y£ 2 \y) + 0(e 3 ). 

This shows that our procedure of averaging-correction is coherent and rigorously justified up to 
order two, for e sufficiently small. 

Extending this kind of reasoning to order greater than 2, following, for example, the procedure 
described in [30, §3], we can show that our procedure is rigorously justified to any order, for 
e sufficiently small, making thereby near-identity transformations. This permits us to obtain 
convergence theorems on n^-order averaging, with suitable assumptions, such as those used in 
[32] . for instance. These kind of results increase precision but not the time-scale of validity. The 
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latter can be improved if the system shows some attracting properties, as it is classically done 
[291 EQl [33]. 

From a practical point of view, we obtain more manageable formulae than those usually 
obtained [10, 32j. Indeed, our approach based on Proposition 13.11 gives an efficient procedure of 
construction of any higher-order averaged systems, a procedure which can be implemented on 
a computer using a symbolic computational software, as it is done for Lie transforms in general 
[21 [251 E] • Actually, many authors [3Ql [31] pointed out the fact that Lie transforms can be used 
to obtain higher-order averaged systems, but in our knowledge, a systematic explicit exposition 
like the one given here has not been published (although a related paper [36J presents some of 
the elements.) The fact that the classical formulae are not manageable may be more evident in 
|15j , where the authors develop an algorithm of higher order averaging for linear equations with 
periodic coefficients, in order to simplify the one of |10| . 

Finally, it is interesting to note that usually the i^-correctors are thought to be the anti- 
derivative in time of the oscillatory part associated to each order (cf. (6) via (7) and (5) of 

|32j for example). This structural form is the same for correctors wffl in our computations, 
u\ 

with Y given by (13 . 3|) instead of the classical averaged part of order i, as a consequence of 
Lie's equations (the reader will be convinced by comparing (6) of |32] for i = 2 and (I3.16p . for 
example) . 

4. Choice of appropriate initial data and the m th approximation of the n th type 

Let n be a positive integer. Consider G P5°(f2), which satisfies the hypothesis of Propo- 
sition 13.11 Then, by E j2.ll and Proposition 12. 1\ every solution x of Y through xq G f2 is in 
P^j-correspondence by (ir o oI t ) teR+ to a unique solution x^ of the n th averaged system 
given by (|3. 12[) in Proposition 13.11 We recall that the latter means 

(4.1) x(t,x ) = (7ro0 lit oX t )(x( n )(t,(vro0-JoZo)(xo))), for alH G R+ 

In a number of applications xq is given, and thus, in order to compute an approximation of x 
based on x^ , we have to determine the initial condition appropriate for "P^^-correspondence, 
that is (-7T o <j)~Q ol )(xo) as shown by (|4.1|) . 

This problem of determination of the inverse transformation for a given Lie transformation 
is not new, and several algorithms have been proposed (see e.g. [HI [35]). We present here an 
approach based on [IE] . 

In order to be consistent with our notations used in (12.12D and (I2.13D . we consider K (resp. 
H ) the generator of the inverse (resp. direct) transformation given by its formal T-series 

(4.2) K(t,i) = T—K n (t,0, (resp. H(t,£) = V -r H n (t,0). 

n>0 n>0 

In the context of our framework, Proposition 5 of [16] can be reformulated as the following 

Proposition 4.1. Let cpi^ generated by Ht defined in \4-<ty an d ^ <fiit generated by the t- 
suspended vector field K t , then: 

(4.3) K t = -{<h,tTH t . 

Proof. We provide a proof for the sake of completeness. Let us denote the inverse of cj> Tt t by rj T) t, 
i.e., for all (r, t) G M x M + , and x G f2, 4>T,t(ilT,t(x)) = x. 
Taking the derivative with respect to r yields 

{D x (t> Tit ){r] Ti t{x))^j^{x) + ^pt(Vr,t( x )) = °- 
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We express the second term on the left-hand side in terms of the generator Ht, and move it 
to the right-hand side to obtain 

(D x <f> Ttt )(r] Tt t(x))^±(x) = -Ht(M T >VT,t(x))), 
which gives, after multiplying both sides with the inverse of (D x (j) Tit )(r} T! t(x)) 

^(x) = -D x ^ t (x) ■ StifoMx))). 

The latter leads to, 

dl, ,t {x) = -D x (f)~U<j> Ttt (r] T>t (x))) ■ H t {<j>T,t{r} T ,t{x))), 



that is, with Definition 2.2, to, 



^i(x) = -(0- t 1 */7,O(/ /7 ,(,)). 



Therefore, for each r and each t, we deduce that 7y T)t is generated by 
(<fi T t)*Ht, and thus substituting r = 1 the proof is complete. 



□ 



The following proposition, which is an obvious corollary of Theorem 12.11 and Proposition 14.11 
where "7/^(1) is given by Definition 2.8, determines formally the inverse transformation for all 
t > 0, and in particular the one in (14. If) for t = 0. 

Proposition 4.2. Let generated by a time- dependent vector field of the form Ht = (l,Ht) T 
and let Kt = —((f>i t)*Ht, then for every non-negative real t: 

(4.4) 7T o 0-1 o l t = Id RN + T Kt (l) ■ K t . 



This proposition allows us to compute the m th approximation of the solution of (jl.ip . which 

lOO 

diff 



is in V^° f ^.-correspondence to one of a n th averaged system by the following algorithm which 



defines the m th approximation of the n th type: 

(1) . Compute the n th averaged system given by (|3.12p and Proposition 13.11 

(2) . Solve Lie's equations as explained in §3.2 until the order m, giving the m first Wt/s 

terms of the r-series of W. 

(3) . The initial condition x(0) of the exact system being given, truncate to the order m the 

ifo-t rans f° rma ti° n (at time t = 0) taken at x(0), Tr: (1) ' -^0(^(0)), where K is loosely 
speaking the generator of the family of inverse transformations, and take this truncation 
plus x(0) as initial condition for the n th averaged system. 

(4) . Compute numerically the solution of the n th averaged system, and add to x^ n \t) 

the truncation to the order m of the series 7wt(l) ■ Wt(x^ n '(t)) , which gives finally the 
m th approximation of the n th type of the exact solution at time t. 

5. Application to a problem in atmospheric chemistry 

The original problem that motivated this work was to examine models of diurnal forcing in 
atmospheric dynamics and chemistry on long time-scales [BJ. The day-to-night changes in the 
radiative heating and cooling of the planetary boundary layer - the lowest part of the atmosphere 
(1-2 km) - are very large. Still, one is often only interested in the slow, season-to-season or even 
year-to-year changes in the way this lower layer interacts with the underlying surface (land or 
ocean), on the one hand, and the free atmosphere above, on the other [8]. The diurnally averaged 
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Figure 1. The figure depicts forcing induced by Si and Zj. 



model we derived here provides insight into a similar problem, that of the basic chemistry of 
slow changes, from one day or week to the next, of a highly simplified system of photochemically 
active trace gases in the troposphere (i.e., the lower 10 km of the atmosphere). The system of 
two coupled ODEs we consider in this paper governs the concentration of the chemical species 
CO (carbon monoxide) and 03 (ozone) [20J, 



(5.1) Ml 



= Si(t) - Z 1 (t)x 1 x 2 
# = S 2 (t) + Z 1 (t)x 2 x l - (1 + Z l (t))x 2 + (Z 2 (t) + S 2 (t))±, 



where t — > x(t) = (xi(t),x 2 (t)) T = ([CO](t), [03](t)) T . This system belongs to the general class 
of ODEs systems given by ([Lip . 

In this system, the diurnal forcing is through the functions Si and Zi (see Figure Q] below for a 
typical example), which can have rather complicated shapes. The system (|5.ip is a very simple 
model for air pollution in an urban environment. Changes in the chemistry at the day-night 
transitions and those in the emission of CO and 03 in a city during a day lead to Si and Zi 
which are only piece-wise smooth. 

While it is not difficult to numerically integrate the system as it is by standard methods, the 
forcing carried by functions Si and Zi can have high frequencies and therefore require a time 
step too small compared to the total time of simulation, T max (say, several tens of years). This 
assumes that the goal is to obtain sufficiently smooth and therefore realistic numerical solutions. 

We discuss in this section how well the numerical solutions of the averaged systems, at orders 
one and two, approximate those of the original one. An important part of the numerical work for 
this particular example is to carry out the transformations between the solutions of the original 
and averaged systems developed in £j2j This includes the transformation of initial values of §H 
Overall, the numerical simulations indicate that the correction performed by Lie transforms to 

(2) 

the first order of the first averaged system M2 or the second averaged system M2 obtained 
by applying results of $3j provides a good approximation to the original one. Furthermore, 
regularity of the averaged systems in time allows one to integrate them by simple methods, such 
as an Euler or a Runge-Kutta of order 4 (RK4 in the sequel). 



20 



M.D. CHEKROUN, M. GHIL, J.ROUX AND F. VARADI 



Last but not least, we can also use larger step sizes for the averaged systems than for the 
original one, which is the key to speed up simulations of long-term dynamical phenomena. 

5.1. First and second averaged systems analysis. Let := M 2 \{x2 = 0}. Then the vector 
field Y associated to Ml belongs to V 00 ^) and Y G C°°(r2). Assume that there exists, for each 
t > 0, a diffeomorphism <f>\f G V%°(fl), generated by G t = {l,Gt) T , such that (4>ij)*Y = Y ave . 

In order to demonstrate the numerical efficiency of the averaged systems, we present only an 
approximation of the family of diffeomorphisms (4>i,t)teR+ U P t° order one. 

For that, we compute the first corrector G °\ = Gq^ (see Definition 2.8) given by integration 
of (|2.22p . Define Go t i(t,^) (resp. Go :2 (t,^)) as the first (resp. second) component of Go(i,£). 
Then, from (|5.ip . 

t ft 

^ ' ' ' ' 2) 



(5.2) Go ,i(t,0= [ SS 1 (s)ds-( 1 4 2 f <JZi(a)<fe + q,,i(£i,£ 

Jo Jo 

(5.3) G ,2(t,0 = 



£ - 1). £ 5S 2 (s)ds + (a - 1)6- Jo SZ x {s)ds 
+^J*SZ 2 (s)d S + C , 2 (^^ 2 ) 



where, for i G {1,2}, 

(5.4) 5Si(t) = Si(t) - Si and SZ^t) = Z l (t) - ~Z iy for all t G 



The constant parts Cb,i and Cq j2 of (|5.2p and (|5.3p are given by assuming that for all £ G Q, 
So Jo Gq(^, s)dsdt = 0, which allows us to avoid "secular" terms at the first order, namely 



T t T t 



(5.5) 0),i(£i,6) = - J J dS 1 (s)dsdt + ^ 2 J J SZ l {s)dsdt, 



and from (15.3 



(5.6) C , 2 (£i,6) 







<i - l)J J SS 2 (s)dsdt - (6 - 1)6 JJSZ 1 (s)dsdt 



T t 

-f / / 5Z 2 {s)dsdt 





Now, let x(0) be the initial condition of a solution x in system M2. In order to compute 
the first approximation of the first type, following the procedure described at the end of £jU 
we have to truncate to order one the .Ko-transformation at time t = of Kq taken at x(0), 
i.e. Tr: (1) • Kq(x(0)), where K is, loosely speaking, the generator of the family of inverse 
transformations. Thereby, we get a first approximation of the initial condition x(0) of the 
solution x which is 'P^^-correspondent to x by the family (<f>i ; t)teM+ ■ 

So we have by Proposition 14.21 end of §3] and Theorem 7.1, that the first approximation of 
the first type x^^^O) of the initial condition x(0) is 

(5.7) x^(0) = x(0)-G (0,x(0)), 

with Go(0,x(0)) determined by the preceding constants evaluated at x(0). Note that the first 
index of the exponent (1,1) in (15.7P indicates that we deal with the first type of approximation 
and note that the second superscript indicates that we make a truncation up to order one in 
the series Tk (1) ■ Kq(x(0)). An identical convention will be used in the following when we work 
with approximations of higher-order type. 
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Next we numerically compute the solution z of M2 through x^ 1 ' 1 ' ) (0) and then we define the 
first pullback x^ by truncating to "order one" the expression 
z(t) + T Gt (l) • G t (z(t)), which gives: 

(5.8) x (1 \t) : = z(t) + G (t,z(t)), for all t G R+. 

We describe now the construction of the second pullback. By Proposition 13.11 with n = 2, the 

second averaged system M2^ 2 corresponds to the vector field = Y + \Y^ where Y^ is 
determined according to (|3.3|) with m = 2. 

(2) 

Both direct computation by hand and symbolic manipulation software for Y$ 
yield for system (|5.1j) the following expression as function of £i and £ 2 , namely 

' t T t 

{fSS^ds - J J SSi^dsd^.i-b-o-Z^t)) 
00 

t T t 

+(j5S 2 (s)ds -J J 5S 2 (s)dsdt).[Zx.(l - h)-oZx{t)\ 

00 

1 T t 

+(J SZx(s)ds - J J SZ^dsdt) < 
00 



v m v (2) 
2 0,1 ' 2 0,2 



(2) 



0,1 



T 

— / < 







+^f 2 (i) 



■I, 



T t 



+ (/<5Z 2 (s)ds -Jf6Z 2 (s)d8dt).(-&.<TZ 1 {t)) 







f2 



and: 



(2) _ 

0,2 _ rp 



T 

- I < 







(j8Si(s)d8 - J $l5S x dsdt).(£ 2 .o-Z x {t)) 


( [aaZi(t) 

t T t 

+{f5S 2 (s)ds- f j5S 2 dsdt). < 
00 " 



T t 



+^Z 2 (t) 
I +(3aS 2 (t)] 



> dt, 



+(f6Z 1 (s)ds-Jj6Z 1 dsdt). { +i°S 2 {t) 







t T t 

+(J 5Z 2 (s)ds - J / SZ 2 dsdt). < 
00 " 



(fr-i)£aZ 2 (t)] 

--] 



where a (resp. /3, 7) are defined by a = 2^- — — 6 + 1 (resp. /3 = 7£~p"(^ — 1)) 7 = 
(6 ~~ ~~ r"))> °~Zi(t) = Zi(t) + Zi and aSi(t) = Si(t) + Si for all t. The terms SZi and 55j 



6' 

are defined in (|5.4p . 

Denote by x^ 2 ^(0) the initial condition for M2 2 through which the solution x^ is 'P^ff 
correspondent to x. We take as approximation of x^ 2 \0), following the preceding procedure and 
notations concerning the first averaged system and taking into account the invariance of the first 
corrector (cf (|3.18p ). the vector 3?( 2,1 )(0) is given by 

(5.9) ^ 2 ' 1 )(0) = x(0)-G (0,x(0)), 

i.e. ^^(O) defined in (|57j| . 
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comparison between the 1st, 2nd avraged solution & the full 

1.8 
1.6 
1.4 
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0.8 
0.6 



0.4 
0.2 

2 4 6 8 10 12 14 16 18 20 

Figure 2. The figure shows the comparison between the full solution and the 
I s * and the 2 nd averaged solutions, on the 2 nd component, for the total time of 
simulation. 

Then we numerically compute the solution v of M2^ based on x^ 2 '^(0), and define the second 
pullback as the first approximation of the second type, i.e.: 

(5.10) x (2) (t) := v(t) + G (t,v(t)), for all t € R + . 

The choice of computing only the first approximation, based on the first and second averaged 
systems, allows us to compare the approximations of the original system M2 by these averaged 
forms, the corrector being the same in each case, but applied to different solutions. This will be 
discussed in §5.31 

The numerical experiments we present here correspond to choices of the forcing functions 
(Si, Zi) which produce "broad" oscillations in the solutions of the original system. These are not 
very realistic but our goal here is to test how well the method performs for "large perturbations" . 
The magnitude of perturbations is shown in Figure [TJ In all the numerical simulations, Z\ and 
S\ were oscillatory but we used the constants Z2 = 5.0 x 10~ 2 and S2 = 12.0 x 10~ 2 . 

5.2. Gain in numerical efficiency. The original system is integrated by standard methods 
(Euler, Runge-Kutta,...) with step size St, which has to be small enough to resolve the oscilla- 
tions induced by the forcing terms Zi and S{. The averaged system M2 or M2^ is integrated 
with step size At, which is typically ten times larger than 5t. 

A key point in numerical efficiency is the regularity of the solutions of the averaged system. 
Indeed, if the solutions of the averaged system are smooth enough in the sense that they do not 
show multiple oscillations over one forcing period, which is one day in our case, we can integrate 
the averaged system with a large step size At, without loss of regularity. This fact for system 

(|5.ip is pointed out in Figure [2j for instance, where the solutions of the averaged systems M2 
(2) 

and M2 , on the second component, are those which do not show multiple oscillations over one 
day. This fact is well known in the e-dependent case where the drift described by the averaged 
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system can be integrated with a step size chosen to be \ times larger than for the non-averaged 
system [2]. 

In order to investigate how well the solutions of the averaged system corrected up to order 
one approximate that of the original one, we did the following: 

(i) For a given initial condition for the original system, x(0), we compute the l s *-approximation 
of the l si -type ^ 1 ' 1 '(0) and the l st -approximation of the 2 nd -type aK 2,1 )(0) given by (|5.7p 
and (|5.9p . of the initial conditions x(0) and x^ 2 ^(0), respectively. 

(ii) Starting from these initial conditions, we compute x and x^ by integrating M2 and 

(2) 

M2 with At as step size, by using a standard integrator. This costs less than in- 
tegrating M2 with a small enough step size that resolves the oscillating forcing terms. 
Here it is sufficient to use the Euler method. 

(iii) We do a simple linear interpolation on x and x^ to obtain the solution on the temporal 
grid defined by St. 

(iv) According to ()5.8[) (resp. (|5.10p ) and the fact that the correctors are only composed of 
integrals, ()5.2[) and ()5.3[) are used for computing, on the grid defined by St, the first and 
the second pullback based on solutions obtained in (ii). 

The numerical tests indicate that the solution of M2, obtained by the procedure described 
through the items (i) to (iv), has accuracy which is close to that obtained by integrating M2 
itself by classical methods, such as a RK4 scheme. The CPU time linked with our procedure 

depends essentially on the one defined for solving the averaged systems M2 or M2^ with step 
size At, and on the one for the integral correction procedure. 

We observed that the use of the first pullback allows to obtain a gain in CPU time nearly 
equal to 75 percent when we use a large step size for At, while retaining good accuracy and 
regularity of the solutions, as shown in the following subsection where the tests and the related 
figures are for St = 0.01 and At = 0.1. 

5.3. Accuracy achieved by the first and second pullback. One of the main reasons for 

(2) 

computing the second averaged system M2 were to obtain a better approximation. In this 
section we demonstrate numerically this fact for St = 0.01 and At = 0.1. Note that the tests 
performed by RK4 on M2 with At give no smooth solutions (not shown), whereas, as shown in 
Figures [3] and U this is not the case for solutions obtained by our method described in (i)-(iv) 
of 

Based on the method of approximation used here, in many cases, the second pullback is better 
than the first, and gives a solution which is close to the one obtained by RK4. In Figures [3] 
and [U we observe that we have coincidence of the I s *-, 2 nd -pullback and the full solution for 
the global trend in each component. A more accurate analysis shows that the 2 nd -pullback is a 
better approximation than the first as shown in Figure [5] and Figure [6] where the 2 nd -pullback 
is represented by dash, the l s *-pullback by dots, and the solution computed with RK4 by solid 
line. This fact is made more evident by the analysis of the error in the supremum norm as shown 
in Figure where the error produced by the 2 Tld -puh 1 back is under 6 x 10 -3 and in Figure [8] 
where the error with the 2 nrf -pullback is represented by dash. 

Nevertheless, we can easily imagine that for a system of dimension larger than two the CPU 
time required to calculate the 2 nd -pullback would increase with algebraic complexity of the 2 nd - 
averaged system (cf. §5.ip , but as we shall see in g6]the 2 nd -pullback can be used for the problem 
of determination of periodic solutions in a general framework. 

The tests presented in this study give satisfactory results, an accuracy of order 1 x 10 -3 for 
the second pullback. Moreover, we have to note that this accuracy is obtained with oscillations 
actually "far" from their averages (see Figure [Q again) producing oscillations, of magnitude 
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comparison between the 1 stthe 2nd Pullback and the full on the 1 st component 




2 4 6 8 10 12 14 16 18 20 

Time(Days) 



Figure 3. The figure shows the comparison between the full so- 
lution and the I s * and the 2 nd pullback solutions, on the I s * com- 
ponent, for the total time of simulation. 

comparison between the 1st,the 2nd Pullback and the full on the 2nd component 

1.8 i 1 1 1 1 1 1 1 1 1 1 
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Figure 4. The figure shows the comparison between the full so- 
lution and the 1 st and the 2 nd pullback solutions, on the 2 nd com- 
ponent, for the total time of simulation. 
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comparison Detween tne isune ^na fuiiDacK ana tne tun on tne ist component 



2nd pullback 

full solution 

• 1st pullback 




17 17.5 18 18.5 19 19.5 

Time(Days) 

Figure 5. The figure shows the comparison between full solution 
and the I s * and the 2 nd pullback solutions, on the 2 nd component, 
for the periodic regime. 

comparison oetween tne isune ^na rullDacK ana tne tun on tne zna component 




18 18.5 
Time(Days) 



Figure 6. The figure shows the comparison between full solution 
and the I s * and the 2 nd pullback solutions, on the 2 nd component, 
for the periodic regime. 



about 1 x 10 _1 , on the full solution obtained by RK4, with a significant gain in CPU time 
compared to the standard method. 

Therefore we can conclude that the Lie transforms averaging method developed here gives a 

(2) 

rigorous setting for constructing the correctors and M2 to provide numerical approximations 



26 



M.D. CHEKROUN, M. GHIL, J.ROUX AND F. VARADI 



even for relatively large perturbations in time induced by the forcing terms. Furthermore, 
according to the tests performed in this work, the 2 nrf -averaged system analysis seems to be 
more relevant than the first. 



trrors in tne sup norm ror tne 1st component 




error with 2nd pullback 

— error with 1 st pullback 




:*M;>s;i!i iii ill 



10 12 
Time(Days) 



14 16 18 20 



Figure 7. The figure shows the error, in the supremum norm, 
between the full solution and the 1 st and the 2 nd pullback solutions, 
on the I s * component, for the total time of simulation. 
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Figure 8. The figure shows the error, in the supremum norm, 
between the full solution and the 1 st and the 2 nd pullback solutions, 
on the 2 nd component, for the total time of simulation. 
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6. Computing periodic solutions by higher-order averaged systems analysis: 

some considerations 

In the classical approach, i.e. in the e-dependent case, in order to find periodic solutions 
the theory of averaging can be a useful tool as shown, for instance, by theorem 4.1.1 (ii) of 
J. Guckenheimer & P. Holmes |13j for the first averaged system, Hartano and A.H.P. van der 
Burgh in [15j . or A. Buica &; J. Llibre in [3] with relaxed assumptions, for the first, second 
and third averaged systems, based on Brouwer degree theory. All these results are linked to a 
maximal size of perturbation eo under which the existence of a hyperbolic fixed point po of the 
first averaged system gives the existence of a unique periodic solution of the original system, 
revolving around pq. 

As was explained in §3.31 the analysis performed in this paper permits to get such results for 
higher-order averaging analysis based on explicit formulae (Proposition 13.11 for averaging and 
Theorem 12.11 for corrections) by revising the proof of classical results for first and second order 
averaged systems, given in the literature (e.g. [13]). For instance, if the k (k > 1) first averaged 
systems vanish identically, then a proof of existence and uniqueness of a T-periodic solution in 
a e-neighborhood of a hyperbolic fixed point of the (k + l)* /l -averaged system can be given on 
the basis of principles given in the proof of theorem 4.1.1 of [13], essentially by showing that the 
Poincare maps of nonautonomous and autonomous systems are e-closed. 

In the light of such considerations, a question arises naturally: do there exist such results for 
the e-independent case? 

A preliminary analysis can be done on time T-maps. Indeed, we can formulate the following 
lemma: 

Lemma 6.1. Let Y £ T' r (£i), and Z € C' r (f2), (1 < r < oo). Then Y and Z are V^j j -equivalent 
if and only if their time T-maps are C r -conjugate. 

This lemma gives thereby a way to study sufficient conditions for the existence of solutions of 
the nonlinear functional equation (|2.6p . which will be investigated in a forthcoming paper. Note 
that the existence of time T-maps is realized under the assumption (A). 

In practice, if the conjugacy of time T-maps is achieved, we can localize a T-periodic orbit 
of the nonautonomous system Y. Indeed, suppose that the time T-map associated with an 
averaged system Z has a fixed point. Then by conjugacy, this fact still holds for the time T-map 
associated with the nonautonomous system. Making use of the proof of Lemma 6.1, we can 
observe that if (tt o c/> ljt o Z t ) t&R + denotes the solution of (|2.6p obtained by Lie transforms, then 
7ro 0^^0X0 is the diffeomorphism realizing the conjugacy between time T-maps associated with Y 
and Z. Therefore Proposition 22] and Theorem 12.11 allow us to compute an approximation of this 
diffeomorphism, which can be applied to a fixed point rj associated with the system Z, leading 
to an approximation £0 of an initial datum lying on the T-periodic orbit of the nonautonomous 
system which is in T^j-correspondence with rj. Thus, by a standard integrator based at £o> we 
can compute an approximation of this T-periodic solution. Such a method can be a useful tool 
for localizing T-periodic solutions of T-periodic nonautonomous dissipative systems, which are 
known to exist, in this case, as shown by classical results on the topic (cf. [271 P- 235]). 

This procedure is relevant for our system (|5.ip . as proved by elementary analysis. Indeed, 
there exists a hyperbolic fixed point po = (1.584, 0.431) T of M2 and, as we can see on Figures El 
[5] and Figures [U [U the first pullback solution is a good approximation of the periodic solution 
revolving around the fixed point. 

Finally, note that the notion of equivalence considered in this paper (Definition I2.ip is the 
appropriate one from the numerical perspective described here in order to compute T-periodic 
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solutions, in view of the fact that the main ingredient was to achieve a correspondence be- 
tween fixed points (with all the periods) of an autonomous system and T-periodic orbits of a 
nonautonomous one. 

7. Appendix 1: Solution of the time-dependent pullback problem via Lie 

transforms 

Let p be an integer greater than or equal to one. Let A(r, x) be a smooth vector field on 
RxP expanded in powers of r as 

(7.1) A(r,x) = A T (x) = Y,-A^(x), 

^-^ nl 

n>0 

and let H (r, t, x) be a smooth vector field on R x R + x W expanded as, 

(7.2) H (r, t, = H T (t, = H T)t (0 = ]T ^H n (£, t) = £ ^#n,t(0, 

n>0 ' n>0 

where t is fixed. Denote by 3( the semiflow generated by Ht, namely the semiflow of 

(7.3) s Ht { £ = h(tM) ;^»VeK. 

Define the Lie transform generated by Ht of A evaluated at r, denoted by L(i?t)(r) • vl, as the 
vector field: 



(7.4) L(tf t )(T)-A = 4°J + £^4?. 



m>l 

where the sequence of vector fields (acting on MP) {A^}, is calculated from the sequence {^4^} 
given by (|7.ip . using the recursive formula: 

n 

(7.5) Alf ] = A« M + ^ (7^_ fe>t 41; for all (i, n) G Z* , 

fc=0 

with Ln n _ k t expressing the Lie derivative with respect to H n _k,t (see e.g. [35] , for details). 

A classical result is the following theorem which expresses the pullback of a vector field as a 
Lie transform which can be proved easily, by adapting, for instance, the work of [19] relying on 
Taylor series expansion and the key formula 

(7.6) ^r,t)*A T = (Z T)t )*(L Hrt A T + d T A T ) 

UT 

from the realm of differential geometry (e.g., [T| I28j). 

Theorem 7.1. Let r G R, as above. The pullback at time t of A T by H Tj ( generated by Ht is the 
Lie transform generated by Ht of A, evaluated at t, that is: 

(E T:t )*A T = L(Ht){r) ■ A, for all t G R+. 

Let B be another smooth vector field on R x W with formal expansion 

(7.7) B T = V — B m . 

ml 

m>0 

One of the advantages of considering pullback as a Lie transform, is that the framework of the 
latter yields linear conditions for finding E T)t which satisfies (H Ti t)*^4 T = B T . Indeed, using the 
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formal power series expansions (|7.4p and (|7.7| . (H T)4 )*.A r = B T leads to a sequence of recursive 
linear PDEs (e.g. [35]), namely, 

(7.8) A$ =B m ,m£ Z+, t G R+ 

where the (n < m) present in (j7.2j) are the unknowns contained in An^ , determining by 
this way the generator Ht of the Lie transform. 

These equations are usually called Lie's equations. 

As a consequence, we can state the following corollary of Theorem 7.1: 

Corollary 7.2. A necessary condition for the existence of a two-parameter family of diffeomor- 
phisms (^T ) t)( T ,t)e'RxM+) generated by the one-parameter family of vector fields (-Ht)t e a+ given 



by JZi^p , such that (z, T j)*A T = B T for all t > and all rel, is that Lie's equations ( 7.8) are 
solvable for the unknowns H m j for all (m,t) G Z + x M + . 

8. Appendix 2 (added for the ArXiv version): Proof of Lemma 6.1 of Chekroun 

et al, DCDS, 14(4), 2006. 

We give here the proof of the Lemma 6.1 of Chekroun et al, DCDS, 14(4), 2006, that 
corresponds here also to Lemma 6.1 of the present manuscript. 

Proof, of Lemma 6.1. As the vector field are assumed to be complete on 0, there exists an open 
subset Vc!J, such that the time T-map associated with Z, denoted by P is well-defined on V. 
We denote by P the time T-map associated with the periodic vector field Y. 

First of all, we suppose that Y and Z are P^^-equivalent, so there exists a map $ G P|(Q), 
such that: 

(8.1) x(t,x ) = * t («(t ) *Q 1 (a;o))) J for all ,x G O, 

using the notations of Definition 12.11 

Since for all (6 V, P(£) = x(T,£) and P(£) = z{T,£), then we get from (j8Tl) : 

P o $ (x ) = x(T, $ (xo)) = *t(^(T, x )) = $ Q (z(T, x )) = $o o P(x ), 

and the necessary condition of the lemma is satisfied. 

For the sufficient condition, let us denote by zt ■ zq — > z(t, zq) the flow of Z such that 
z(0, Zq) = zq, and by x t : xo — > x(t,xo), the semiflow of Y such that x(0, xo) = xq. Then 
zt+T = ztoP (at least in V) and xt+r = xt°P. By assumption there exists a C r -diffeomorphism 
H such that P o H = H o P. Let us introduce for all i > 0, 

H t = x t oHo( Zt )- 1 . 

Then Ht is a time-dependent C r -diffeomorphism with Hq = H. Obviously Ht carries solution 
z(t,zo) into x(t, H~ l (z$))\ and we have: 

(8.2) H t+T = x t+ T o H o {zt+r)- 1 = x t oPoHop- 1 o {z t )- 1 = x t o H o (z^ 1 = fit, 
that gives the T-periodicity of the change of variables. □ 
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